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The method of wavefront shaping to control optical properties of opaque media is a promising tech¬ 
nique for authentication applications. One of the main challenges of this technique is the sensitivity 
of the wavefront-sample coupling to translation and/or rotation. To better understand how trans¬ 
lation and rotation affect the wavefront-sample coupling we perform experiments in which we first 
optimize reflection from an opaque surface - to obtain an optimal wavefront - and then translate or 
rotate the surface and measure the new reflected intensity pattern. Using the correlation between 
the optimized and translated/rotated patterns we determine how sensitive the wavefront-sample 
coupling is. These experiments are performed for different spatial light modulator bin sizes, beam 
spot sizes, and nanoparticle concentrations. We hnd that all three parameters affect the different 
positional changes, implying that an optimization scheme can be used to maximize the wavefront- 
sample coupling’s stability. We also develop a model to simulate sample translation/rotation and 
its effect on the coupling stability, with the simulation results being qualitatively consistent with 
experiment. 


I. INTRODUCTION 

In 1990 Freund predicted that precise optical de¬ 
vices could be made using opaque media and wave- 
front shaping [l| . Since then this technique has been 
used to control: transmission through opaque mate¬ 
rials , the polarization of light [j, 0 , broadband 
spectral characteristics , and the spatio-spectral 
properties of random lasers |ili3. It has also been 
used to enhance fluorescence microscopy [a 
achieve perfect focusing (ailllj compress ultrashort 
pulses [1 ^, [itI a nd enhance astronomical and biologi¬ 
cal imaging [isL El- Recently the technique’s biolog¬ 
ical applications have been supplemented by photo 
acoustic wavefront shaping (PAWS) which combines 
wavefront shaping and ultrasound technology to bet¬ 
ter image biological tissue pol-l^ . 

Additionally, using optimal wavefront shaping to 
control optical properties of opaque media has been 
proposed as a method of implementing optical phys¬ 
ically unclonable functions (PUFs) |26l-[^. PUFs 
are systems which use randomly distributed proper¬ 
ties of a material (scatterers in opaque media, de¬ 
fects in semiconductor chips, etc) to create unique 
unclonable signatures [^, The unique signa¬ 

tures of PUFs are of particular interest in authenti¬ 
cation applications (such as authenticity verification 
and tamper indication) as their uniqueness makes 
counterfeiting a PUF unfeasibly difficult, if not im¬ 
possible. 

The envisioned application of optimal wavefront 
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shaping to the implementation of optical PUFs is as 
follows: an opaque surface marker - either intrinsic 
to the surface (such as defects or roughness) or an 
extrinsic glue (with scattering particles) - is inter¬ 
rogated using a shaped wavefront and the surface’s 
optical response is measured. The response is then 
optimized using an optimization algorithm in order 
to determine an optimal wavefront. The optimal 
wavefront and surface response are then stored as a 
challenge-response pair, which is used at later times 
to verify that the surface is unchanged. Changes in 
the challenge-response pair would indicate that the 
surface had been tampered with. 

In an ideal implementation of this type of op¬ 
tical PUF, the changes in the challenge response 
pair would occur only due to intentional tampering. 
However, in reality there are other unintentional ef¬ 
fects which can change the challenge-response pair 
such as: dust accumulation, material changes due to 
the environment (temperature, humidity, radiation), 
drifts in the optical equipment, and sample position¬ 
ing changes. These various effects can work together 
to decouple the challenge and response, leading to 
a false indication of intentional tampering. Given 
these concerns it is necessary to quantify these un¬ 
intentional effects if this technique is to be used in 
the field. 

In this study we specifically quantify the effects 
on the challenge-response coupling due to changes 
in sample positioning, including translation and ro¬ 
tation. For our PUF samples we use a nanoparticle 
(NP) doped polymer and for the optical response 
we measure the reflected intensity pattern, with op¬ 
timization working to focus the scattered light into 
a target area. In addition to experimentally measur- 
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ing the effects due changes in sample positioning, we 
also use an extended random phase Gaussian beam 
model [13 to simulate translation and rotation in 
order to better understand the underlying physical 
mechanisms. 


II. THEORY 

A. Extending the Random Phase Gaussian 
Beam Model 

Previously, we developed a random phase Gaus¬ 
sian beam model (RPGBM) [13 to model optimized 
transmission through an opaque material. One of 
the weaknesses of this model was the treatment of 
the material as only affecting the phase of the inci¬ 
dent light, and not the amplitude. In order to take 
these effects into account we use a random matrix 
theory (RMT) model of scattering in disordered me¬ 
dia [3, IsiMioj . In the RMT model, scattering in a 
disordered media is represented by complex-valued 
matrices which map the incident electric field at in¬ 
put channel m to the exiting electric field at output 
channel n. Specifically, using the RMT formulation 
of Vellekoop and Mosk, we can write the electric field 
at the output channel n as 


N 


-^t;n — ^ ^ tnmEm’) 

(1) 
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N 
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IM 
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with Em being the incident electric field at input 
channel m, are the transmission matrix ele¬ 

ments, which maps the input channel to the 
transmitted channel, and rnm are the reflection 
matrix elements, which map the input channel 
to the reflected channel. In this study we define 
the input and output channels to be related to the 
spatial coordinate on the sample surface, with the 
coordinates discretized for calculation purposes. 

In Vellekoop and Mosk’s model they only consid¬ 
ered the transmission matrix elements [40| , while in 
this study we are concerned with reflection. There¬ 
fore at this point we will derive the reflection matrix 
elements. Assuming no absorption, the transmission 
matrix, t, and reflection matrix, r are dependent on 
each other via 

rtr + t1't = I, (3) 

where I is the identity matrix and f denotes the con¬ 
jugate transpose. The form of Equation [3] arises due 


to conservation of energy in the absence of absorp¬ 
tion. From Vellekoop and Mosk’s model we know 
that the transmission matrix can be modeled as a 
random unitary matrix, U, drawn from a circular 
Gaussian distribution and scaled by a scalar trans¬ 
mission coefficient, t, giving the transmission matrix 

Hill: 

t = tv, (4) 

where the form of Equation [4| is a consequence of 
treating the scattering system using RMT. 

Substituting Equation [4| into Equation [3] and re¬ 
arranging terms we can derive an equation for the 
reflection matrix 

r1’r = I-I^I^WU, (5) 

rtr= (1-|t|2)I, (6) 

which can be solved with a reflection matrix defined 
as 


r = rV, (7) 

where V is a unitary matrix and r is a scalar reflec¬ 
tion term satisfying, 

= (8) 

Since the reflection matrix terms should be indepen¬ 
dent for a disordered medium [32|, [H, [H, |i3 , for our 
modeling we can draw the random unitary matrix V 
from a circular Gaussian distribution. This choice is 
also consistent with the maximal fluctuation theo¬ 
rem, which predicts that all eigenchannels are either 
open or closed [H, HI, [13 • 


B. Simulations 

To model the effect of translation and rotation on 
the coupling between an optimal wavefront and sam¬ 
ple we begin by assuming a Gaussian beam incident 
on the sample with an electric field of: 

Em = Eo exp |_ (^ ^yo) -^ 

where mo is the peak position, a is the Gaussian 
width, and ^rn is the phase mask due to the SLM. 
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Substituting Equation [9] into Equation [2] and apply¬ 
ing the RPGBM we find the electric field in the de¬ 
tector plane to be: 

N N 

Ed-n' = EE rnmEo exp 

n=0 m=0 

-i^nn' - iV’m 

( 10 ) 

where Ax is the computational grid spacing, N is the 
total number of grid spaces, and the term i^nn' is 
the exponential term in the Discrete Eourier Trans¬ 
form. Using Equation [TO] and a microgenetic algo¬ 
rithm [HI we systematically vary ^rn and calculate 
the intensity in the detector plane 


(m — mo)^Ax^ 



m = 1 


m = 16 


EIG. 1: (Golor online) Schematic diagram of 
transverse translation modeling using a grid with 
16 points and a bin size of 6 = 2. As the 
translation. Am, is increased the wavefront shifts 
to the right, with the left wavefront bins becoming 
zero, and the right wavefront bins going off grid. 
Note that while not shown the amplitude shifts as 
well. 


— \Ed-,r> 


( 11 ) 


C. Translational Stability 


to find the phase mask giving a focused spot in the 
detector plane. At this point we note that during 
experiments we typically bin the SLM pixels into 
super-pixels of side length b (i.e. the bin area isbxb 
px^). Eor our 1-D simulations we still use bins but 
with them corresponding to groupings of grid points. 

Once the optimized wavefront is found we use nu¬ 
merical transformations, discussed below, to model 
sample translation and rotation. Eor each positional 
change we recalculate the intensity profile in the de¬ 
tector plane and compare the new profile, J^, to the 
original optimized profile, /o,i using the Pearson cor¬ 
relation coefficient: 


c = 




i=0 


N _ N 

EX-irEiK 

i=0 i=0 



( 12 ) 


where /o,i is the initial optimized profile, U is the 
intensity profile after translation or rotation, /q is 
the average value of /o,i, and / is the average value 
of li. The coefficient varies from -1 to 1 with 1 cor¬ 
responding to perfect correlation and -1 correspond¬ 
ing to perfect negative correlation. When perform¬ 
ing simulations we use one spatial dimension, due to 
computational restraints, and Equation [12] is used 
as written with only one index. However, in our 
experimental results we use two dimensions, which 
requires Equation [12] to be extended with a second 
index and all sums are turned into double sums. 


1. Transverse Translation 

Once the optimal wavefront, has been deter¬ 
mined and the optimal intensity profile recorded we 
can model transverse sample translation by shifting 
the incident field’s index by Am, changing the elec¬ 
tric field in the detector plane to be: 

N N 

Ed\n' = EE ^nmEo exp 

n=0 m=0 

(13)' 

The new intensity profile is then calculated using 
Equation [TT] and is compared to the untranslated 
profile using Equation [12] Eigure [T] shows a sketch 
of the wavefront for different translation values with 
the sample grid marked out and a bin size of 2Ax. 

Using this method we calculate the correlation co¬ 
efficient as a function of transverse translation for 
three different bin sizes, shown in Eigure [2l and four 
different incident spot sizes, shown in Eigure [3] The 
model results, as a function of translation, are found 
to follow a stretched exponential. 


(m + Am — mo)^ Ax^ 


(7^ 


c(x) = A exp 



(14) 


where cq is a small offset from zero, A = 1 — cq , and 
the functions width, ax, and exponent, /3 depend 
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Translation (units of Ax) 


FIG. 2: (Color online) Simulated correlation 
coefficient as a function of sample translation for a 
Gaussian width of 200 Ax and three different bin 
sizes with fits to a stretched exponential function. 
As the bin size decreases the width of the stretched 
exponential decreases. 


on the bin and spot size. Equation [14] has a Half- 
Width-Half-Max (HWHM) of W = cr^(In 

From Figure |2] we see that as the bin size increases 
the HWHM of the correlation curve increases imply¬ 
ing greater resistance to decorrelation as the sample 
is translated. This can be understood as the frac¬ 
tional translation is inversely proportional to the bin 
size (i.e. a translation of 2 Ax is a complete bin 
shift for 6 = 2, while it is only 1/5 of a bin shift for 
b = 10), meaning the change due to translation is a 
smaller effect on the wavefront. 

The other factor determining the fractional trans¬ 
lation amount is the spot size of the incident beam, 
with the fractional translation being inversely pro¬ 
portional to the beam width. Figure [3] shows the 
correlation curves for different spot sizes, with the 
larger spot sizes resulting in greater HWHMs. As 
with bin size, we find that the larger the segments 
of the phase front, the more resistant to transverse 
translation the wavefront coupling becomes. While 
these results suggest greater stability for larger bin 
sizes and spot sizes, the benefit needs to be weighed 
against the smaller enhancements that result when 
increasing the bin size or spot size • 



Translation (units of Ax) 


FIG. 3: (Color online) Simulated correlation 
coefficient as a function of sample translation for a 
bin size of 6 = 5 and four different beam spot sizes 
with stretched exponential fits. As the spot size 
decreases the width of the stretched exponential 
decreases, suggesting that smaller spot sizes result 
in less coupling stability. 


2. Translation Along Optical Axis 

In addition to translating the sample in the trans¬ 
verse axis after optimization, we can also translate it 
along the optical axis resulting in the incident spot 
size changing. We model this effect by transform¬ 
ing the beam width, wq and bin size, b gb 

-such that the ratio w/b is unchanged - as well as in¬ 
troducing the effect of the changing Gaussian beam 
curvature given by 


Rig) = zrVg'^ "^(^^ + ^ 2 / 1 ^ ( 1 ^) 

where zr = is the Rayleigh range and k = 

27r/A, with A being the wavelength of the light used. 
Figure m shows a sketch of the wavefront for different 
spot sizes, with the bin size increasing with larger 
spot sizes. Note that Figure [4] does not include the 
phase related to the Gaussian curvature, nor does it 
show the amplitude profile, which also changes with 
spot size. 

Since we experimentally measure the change in 2 ; 
position, not the change in spot size, we can convert 
from the spot size increase, g, to change in 2 ; position 
by considering the Gaussian beam width. 
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w = 20Ax 



m = 1 m = 16 

FIG. 4: (Color online) Schematic diagram of 
modeling translation along the optical axis using a 
grid with 16 points. As the sample is translated 
along the optical axis the on sample beam width 
changes, which changes the coupling between the 
wavefront and sample grid. Namely, as the beam 
width increases, the absolute bin size increases. 


l+(i), (16) 
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Z Translation (units of z^) 

FIG. 5: (Color online) Simulated correlation 
coefficient as a function of sample translation for a 
beam with a Gaussian width of 100 Ax and three 
different bin sizes with stretched exponential fits. 
The correlation curves HWHM increases with large 
bin sizes. 


w{z) = Wo] 


which for a new width, gwo, can be rearranged to 
give the 2 ; position: 


2 : = ZR^/g^ - 1, (17) 

Using the changes in the incident electric field de¬ 
scribed above we model the effect of translations 
along the optical axis on the wavefront coupling. 
Figure [5] shows the correlation coefficient as a func¬ 
tion of 2 : translation for three different bin sizes, with 
the HWHM increasing with larger bin sizes. This 
result is similar to the transverse translation case, 
with the explanation also being that larger bin sizes 
reduce the fractional wavefront change due to trans¬ 
lation. 

We also simulate the effect of spot size on the 
correlation coefficient as a function of 2 : translation, 
with the results shown in Figure [6l Note cr = 10 
Ax corresponds to the beam waist. From Figure 
[6] the spot size appears to have little, if any, effect 
on the correlation curves. This result is due to the 
wavefront-sample coupling depending more on the 
fractional bin size, rather than the absolute size of 
the bin on the sample. In this case, since the bin’s 
fractional size is unchanged by the spot size chang¬ 
ing, the change in spot size shouldn’t affect the 2 : 
translational stability. 


D. Rotational Stability 

To simulate the correlation coefficient as a func¬ 
tion of rotation we begin first by considering how 
rotation affects the reflected laser beam. Due to the 
nature of the transformations required, in this sec¬ 
tion we use a continuous coordinate basis to derive 
the transformations with the understanding that cal¬ 
culations are performed using discrete Fourier Trans¬ 
forms. In this case the reflection matrix, r, becomes 
a response function lZ{x — y) where the reflected field 
is given by 


E'{x) = / dyn{x - y)E{y) (18) 

J — oo 

where E{y) is the incident field. 

To determine the effects of sample rotation we 
simulate rotation of the sample around the y axis 
by an angle 0 as shown in Figure 0 Due to the rota¬ 
tion the beam and sample will intersect at different 2 : 
positions as you move from the center. Namely, the 
difference in the intersection position is given by: 

A 2 :(x) = xtan6>, (19) 

where x is the transverse position in the lab (un¬ 
rotated) coordinates. This transverse position de- 









Correlation Coefficient 


6 



Z Translation (units of Zp) 


FIG. 6: (Color online) Simulated correlation 
coefficient as a function of sample translation for a 
bin size of 6 = 1 Ax and three different spot sizes 
with stretched exponential fits. The correlation 
curves are unaffected by the change in spot size. 


z 
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FIG. 7: (Color online) Coordinates of lab frame 
and sample frame (denoted by *). 


pendent offset will result in the addition of a linear 
phase term to the incident field given by 


A 2 p{x) = kAz{x). (20) 

The offset also affects the Gaussian beam parame¬ 
ters of the incident field giving a transverse position 
dependent Gaussian width of 


w{x) = Wq 


1 + tan^ 6>, 


( 21 ) 


and a transverse position dependent curvature of 


R{x) = xtan6> 


1 + 


x‘^ tan^ 0 


( 22 ) 


Finally, rotating the sample results in coordinate 
transformation for the sample surface. The trans¬ 
formation from the sample coordinates (denoted by 
a superscript *) to the lab coordinates is given by 


X = X* cos6>, 


(23) 


where x* is the transverse position on the sample 
and X is the transverse position in the lab coordi¬ 
nates. Note that Equation [23] implies that for a ro¬ 
tation of 6> = 7r/2, all the sample coordinates trans¬ 
form into X = 0, which is expected as the sample is 
parallel to the optical axis in this configuration. 

With the changes given by Equations [2QH23I we 
can write the field reflected from the rotated sample 


as: 
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E'{x*) = J ds^'R{x^ — s^)Eq exp | 


cos^ 0 


'^ 0(1 + 6>tan^ 0 ) 

cos^ 0 


5 * cos ^ tan ^ 


1 I_ ZR _ 

' cos^ 9 tan^ 0 


— i'ip{s'^ cosO) 

— iks* sin6>|. 


(24) 


At this point applying Fraunhoffer diffraction the¬ 
ory will predict the field in a plane parallel to the 
sample. However, the detector is setup to be par¬ 
allel to the unrotated sample. Therefore to account 
for the coordinate transformation between the sam¬ 
ple and lab coordinates we apply angular spectral 
decomposition to calculate the diffracted field in the 
detector plane 

Angular spectral decomposition begins by first 
transforming the electric field reflected from the 
sample into its spectral representation by taking the 
Fourier transform of Equation [Ml 

/ oo 

(25) 

-OO 

where is the angular frequency in the sample coor¬ 
dinates. Using a spectral coordinate transform 
we can rotate the reflected field into the lab coordi¬ 
nates as 


/ oo 

dxE'{x)e^i^^', (29) 

-OO 

where Z is the distance between the sample and de¬ 
tector. 

Using the model of sample rotation described 
above - with the equations converted to a discrete 
basis for computation - we calculate the wavefront 
coupling stability as a function of sample rotation 
for different bin sizes, shown in Figure |8l and dif¬ 
ferent spot sizes, shown in Figure |9l From Figure 
[8] we find that the correlation curve as a function 
of rotation follows a stretched exponential function 
with the bin size having no effect on the HWHM. 
This effect arises due to the relative magnitudes of 
the changes to the beam due to rotation. For in¬ 
stance, if we assume small angles, such that ~ 0, 
Equation [24] becomes. 


£'(v) = £'cos 0{X ^ sin6>) (26) 

where v is the angular frequency in the lab coordi¬ 
nates. Next we calculate the field in the lab’s spa¬ 
tial coordinates, which requires taking the Inverse 
Eourier Transform of the field and the determinate 
of the coordinate transform’s Jacobian, | J(7i*)| [i^ . 
The Jacobian term is required to preserve energy 
conservation as the coordinate transform is nonlin¬ 
ear [13 . These considerations result in the reflected 
field in the lab coordinates being, 

/ oo 

du*£'{u* cose + [A“2 - u* 2]V2 sin6») 

-OO 

x|J(w*)|e-*2’^“*"^, 

(27) 

where the determinate of the Jacobian is given by 

0, 

77 * 

J{u*) = COS0 --J====sm0. (28) 

With E'{x) calculated we now can calculate the elec¬ 
tric field in the detector plane straightforwardly as 


E'{x^)^ / ds^n{x^ - s^)Eo 

J — oo 

X exp I — 

(30) 

where the only remaining influence of rotation is the 
linear phase term ks'^O. As this term only depends 
on the on sample position, and not the bin size it 
is easy to see why the bin size does not affect the 
rotational stability. To justify this for larger angles, 
>> 0, we note that at these larger angles the 
wavefront-sample coupling is already broken by the 
linear phase term. 

While the coupling’s rotational stability is found 
to be independent of bin size, from Eigure [9] we 
find that the spot size has a drastic effect on the 
coupling’s stability, namely, that as the spot size 
increases the coupling stability decreases. This ef¬ 
fect is a direct result of Equation [191 which shows 
that the scale of the z offset is linearly related to 
the transverse position, meaning that for larger spot 
sizes there will be larger z offsets, resulting in greater 
phase changes. 
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FIG. 8: (Color online) Simulated correlation 
coefficient as a function of sample rotation for a 
Gaussian width of 200 Ax and four different bin 
sizes. All four curves are found to be within error 
of each other suggesting that the bin size does not 
affect rotational stability. 



FIG. 9: (Color online) Simulated correlation 
coefficient as a function of sample rotation for a bin 
size of 6 = 5 Ax and four different spot sizes with 
stretched exponential fits. The correlation 
coefficient’s HWHM is found to increase with 
smaller spot size, suggesting that smaller spot sizes 
are more resistant to sample rotation. 


III. METHOD 


To experimentally determine the coupling sta¬ 
bility of an optimized wavefront under translation 
and rotation we use an optimized reflection setup 
and Zr02 doped polyurethane (PU) nanocompos¬ 
ites. The preparation of the Zr02/PU nanocom¬ 
posites is as follows. Spherical zirconia nanoparti¬ 
cles, Zr02 NPs, are synthesized by forced hydrolysis 
followed by calcination at temperature and time of 
600°C and 1 h, respectively [46|. The Zr02 NPs 
are then hydrophobized by boiling the dispersion 
of the Zr02 NPs (2 mg/mL) in 1 vol% solution 
of n-octadecyltriethoxysilane (Gelest, SIO6642.0) in 
toluene m. The hydrophobized Zr02 NPs are 
dispersed in 10 wt% bisphenol A diglycidyl ether, 
BADGE (EPONTM resin 828, Miller-Stephenson 
Chemical Company, Inc.) solution in toluene at 
a 20 mg/mL concentration with the aid of sonica- 
tion. Toluene acts as a thinner for the otherwise 
viscous BADGE fluid. The remaining BADGE is 
then added to yield the desired weight fraction of 
Zr02 NPs. The mixture is further sonicated until 
a homogeneous mixture is achieved and the toluene 
is evaporated in vacuum. Once the BADGE mix¬ 
ture is prepared an equivalent amount of diethylene 
triamine, BETA (Epikure 3223 curing agent, Miller- 
Stephenson Chemical Company, Inc.) is added and 
the mixture is mixed thoroughly and then poured 
onto a rectangular glass slide (l”xl.5”). The poly¬ 
mer is then polymerized in an 80 °C oven for 2 hours 
with the resulting films having thickness of 1-3 mm 
and a scattering length of approximately 4 //m. This 
implies that multiple scattering events occur within 
the sample when illuminated. 

With the samples prepared we place them inside of 
an optimized reflection setup consisting of: a Coher¬ 
ent Verdi Nd:YV 04 GW laser, a Boulder Nonlinear 
Sciences LCOS-SLM, a Thorlabs CMOS camera and 
various optics. See EignreflQlfor a schematic of the 
experimental setup. The Verdi is operated at full 
power (10 W), where the laser is most stable. Since 
the 10 W output is well above the damage threshold 
of the SLM we use a 90:10 beamsplitter to pick off 1 
W of the laser, which is then expanded and passed 
through a half-waveplate (HWP) polarizer pair to 
further control the beam intensity, as well as main¬ 
tain the correct polarization for the SLM. The beam 
is then reflected off the SLM and focused onto the 
sample using a 20x high working distance objective. 
The resulting speckle pattern reflected from the sam¬ 
ple is then collected by the same objective and re¬ 
flected onto the Thorlabs camera. The camera and 
SLM are interfaced to a computer which uses a mi- 
crogenetic algorithm to optimize the wavefront 
such that the reflected light produces a spot focus 
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FIG. 10: (Color online) Schematic of optimized 
reflection setup. 
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FIG. 11: (Color online) Measured correlation 
coefficient as a function of sample translation in 
the transverse direction for a 25 wt% sample, spot 
size of 2.2 jam and three different bin sizes with 
stretched exponential fits. Like the simulated 
results the fit width is found to decrease as the bin 
size decreases. 


on the camera array. 

To translate the samples accurately we use a Thor- 
labs flexure stage with differential micrometers giv¬ 
ing a translational sensitivity of 500 nm. In order 
to minimize the effect of translational cross talk we 
set the sample zero position to be near the middle of 
the full micrometer range, as per the manufacturers 
instructions. For rotating the sample we use a New¬ 
port precision rotation stage with a sensitivity of 30 
arc seconds (145.4 /irad) 

The procedure for measuring the positional sta¬ 
bility of the wavefront coupling is as follows. First, 
the optimal wavefront, which produces the largest 
intensity in the target area, is determined using a 
microgenetic optimization algorithm [42|. Once the 
optimal wavefront is found the target intensity pro¬ 
file is imaged four times and averaged to find the op¬ 
timized intensity profile. The sample is then trans¬ 
lated or rotated a fixed distance and the new inten¬ 
sity profile is measured four times and averaged to 
find the intensity profile corresponding to the new 
position. This process is then repeated until the 
intensity in the target area is the same as the back¬ 
ground, which corresponds to full decoupling of the 
wavefront. Once the intensity profiles are measured 
they are each compared to the optimized intensity 
profile using Pearson’s correlation coefficient. Equa¬ 
tion [121 


IV. RESULTS AND DISCUSSION 
A. Bin Size 

The first experimental parameter we vary when 
measuring the translational and rotational wave- 
front coupling stability is the SLM bin size. For 
these measurements we use a sample with a NP 
concentration of 25 wt%, three different bin sizes, 
b = {8,12,18}, and a spot size of 2wo = 2.2 fim. 

We begin by measuring the correlation coefficient 
as a function transverse translation for the three dif¬ 
ferent bin sizes as shown in Figure [TTl with the data 
fit to Equation [TTl Erom Eignre [TTl we see that as 
the bin size increases the HWHM of the correlation 
curves increases, as tabulated in Table [H Addition¬ 
ally from EigureUllwe see that the correlation curves 
change shape as the bin size increases with the small 
bin sizes having much steeper decays, correspond¬ 
ing to yd < 1, and the larger bin sizes having more 
gradual decays, corresponding to yd > 1. Both these 
results, the change in HWHM and yd with bin size, 
are consistent with simulation. 

After translating the sample in the transverse di¬ 
rection we next translate the sample along the op¬ 
tical axis. Eigure [121 shows the correlation coeffi¬ 
cient as a function of 2 : position for the different bin 
sizes tested. We find that as the bin size increases 
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FIG. 12: (Color online) Measured correlation 
coefficient as a function of sample translation along 
the optical axis for a 25 wt% sample, spot size of 
2.2 jam and three different bin sizes with fits to 
stretched exponentials. 


the coupling stability increases, with the HWHM be¬ 
coming larger as listed in Table [H The curves are also 
found to intersect near zero as the z translation ap¬ 
proaches a value equal to the Rayleigh range. These 
results, the HWHM increase and curve intersection, 
are predicted by simulations. However, the inter¬ 
section point in the simulations is closer to 2 : = 1.8 
zr, while the experimental intersection is closer to 

^ = Zr. 

At this point it is important to note the relative 
stability lengths of transverse and optical axis trans¬ 
lation. From Table [I] we find that for the smallest bin 
size the transverse HWHM is only 37.3 nm while the 
z HWHM is 1364 nm, corresponding to a scale differ¬ 
ence of 36.6 X. This scale difference, also observed in 
the simulations, implies that the wavefront-sample 
coupling is more sensitive to transverse translation 
than to optical axis translation. The underlying 
mechanism of this difference is related to the dif¬ 
ferent effects each translation produces. Transverse 
translation results in the center of the beam chang¬ 
ing position on the sample surface, thereby effec¬ 
tively shifting the fields index, but leaving the re¬ 
flection matrix’s indices the same. This shift quickly 
breaks the coupling between the optimal wavefront 
and the sample. On the other hand, for translation 
along the optical axis the center of the beam is un¬ 
changed and only the beam width shifts. This means 



FIG. 13: (Golor online) Measured correlation 
coefficient as a function of sample rotation for a 25 
wt% sample, spot size of 2.2 fim and three different 
bin sizes. The correlation curves for each bin size 
are found to be within experimental uncertainty of 
each other, suggesting that bin size does not play a 
role in the rotational stability. 


that the field’s positional indices and the reflection 
matrix’s indices match for a longer translation, re¬ 
sulting in the increased positional stability. 

Next we measure the correlation coefficient dur¬ 
ing sample rotation. Figure [l3] shows the correla¬ 
tion coefficient as a function of rotation angle for 
the three different bin sizes, with all three correla¬ 
tion curves found to be within uncertainty of each 
other. Fitting the curves to Equation [14] we deter¬ 
mine their HWHMs, reported in Table HI and find 
that the HWHMs are unchanged as a function of 
bin size. These results are consistent with simula¬ 
tion where the correlation coefficient as a function 
of angle follows a Gaussian with the width being 
unaffected by bin size. 

B. Spot Size 

The next experimental parameter (that can affect 
the correlation as a function of translation and ro¬ 
tation) which we vary is the beam spot size. These 
measurements are performed using a 25 wt% NP 
concentration sample, a bin size of 6 = 16, and four 
different spot sizes having Gaussian widths of a = 
1.1 /im, 60 //m, 240 /rm, and 540 jam. Note that the 
beam waist for the system corresponds to a = 1.1 
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Bin Size 

X HWHM (nm) 

Z HWHM (nm) 

Rot. HWHM (;urad) 

18 

693.6 ± 9.0 

2844 ± 26 

254 ± 12 

12 

137.8 ± 5.9 

2228 ± 32 

261 ± 12 

8 

37.3 ± 3.9 

1364 ± 22 

265.1 ± 6.2 


TABLE I: Translational and rotational HWHMs for different bin sizes. For both x and 2 ; translation 
smaller bin sizes yield larger HWHM, while the rotation HWHM is unchanged, within experimental 

uncertainty, for each bin size. 



FIG. 14: (Color online) Measured correlation 
coefficient as a function of sample translation in 
the transverse direction for a 25 wt% sample, bin 
size of b = 16, and four different spot sizes. 


yum, or a spot size of 2wo = 2.2 jxm. 

For each spot size we first translate the sample 
in the transverse direction and measure the correla¬ 
tion coefficient at each position, as shown in Figure 
M From Figure [14] we find that as the spot size in¬ 
creases the correlation curves broaden (the HWHM 
are listed in Table [III) 5 which is consistent with simu¬ 
lation. The underlying mechanism of this increase in 
stability is due to the larger spot size causing trans¬ 
verse translations to have a smaller effect on the inci¬ 
dent wavefront, as the fractional change in wavefront 
with transverse translation is inversely proportional 
to the spot size. 

We next measure the effect of changing the spot 
size on the correlation coefficient as a function of 
optical axis position. Figure [15] shows the correla¬ 
tion coefficient as a function of 2 : translation for the 
four different spot sizes, with all four curves over¬ 
lapping within experimental uncertainty. This over¬ 
lap results in the HWHM of each curve, see Table 
im being within uncertainty of each other, which is 
consistent with simulation. As discussed in the sim- 


Scaled Translation (units of z,,) 

0.0 0.2 0.4 0.6 0.8 1.0 



FIG. 15: (Golor online) Measured correlation 
coefficient as a function of sample translation along 
the optical axis for a 25 wt% sample, bin size of 
6 = 16, and four different spot sizes. 


ulation section, this effect is due to the bin size to 
spot size ratio being unchanged with changing spot 
size. This implies that the more important quantity 
to the translational stability is the fractional bin size, 
and not the absolute bin size. 

After measuring the translational coupling stabil¬ 
ity with different spot sizes we move on to measuring 
the effect of the spot size on the correlation coeffi¬ 
cient as a function of rotation as shown in Figure [161 
From Figure [16] we find that as the spot size increases 
the correlation curves narrow, see Table El for the 
HWHMs, and the initial decay steepens, with the 
stretch parameter approaching (3 = 1. These results 
are consistent with simulation, with the underlying 
mechanism once again related to larger spot sizes 
producing larger wavefront distortions than small 
spot sizes. 
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Spot Size (/im) 

X HWHM (nm) 

Z HWHM (nm) 

Rot. HWHM (;urad) 

1.1 

270.6 ± 9.4 

2045 ± 40 

462 ± 10 

60 

557 ± 23 

2048 ± 50 

411 ± 12 

240 

699 ± 81 

2040 ± 25 

219.8 ± 7.3 

540 

796 ± 28 

2050 ± 44 

140.9 ± 3.5 


TABLE II: Translational and rotational HWHMs for different incident spot sizes. 



FIG. 16: (Color online) Measured correlation 
coefficient as a function of sample rotation for a 25 
wt% sample, bin size of 5 = 16, and four different 
spot sizes. 


C. Nanoparticle Concentration 

Thus far we have considered the effects of bin size 
and spot size on the wavefront-sample coupling’s 
ability to withstand sample rotation and transla¬ 
tion, with the results found to be consistent with 
simulations using the extended RPGBM. The final 
experimental parameter we consider is the NP con¬ 
centration of the samples, which is not modeled by 
the extended RPGBM. Therefore we do not have 
any simulation results to make comparisons too. In 
order to effectively model the effect of NP concen¬ 
tration the transmission and reflection matrices need 
to be calculated in such a way as to take scatterer 
concentration into account, which cannot be done 
by the random unitary matrix model used by the 
extended RPGBM. Work is currently underway to 
address this limitation, but is beyond the scope of 
this study. 

To measure the effect of NP concentration on the 
wavefront-sample coupling’s rotational and transla¬ 
tional stability we use three different samples with 



FIG. 17: (Golor online) Measured correlation 
coefficient as a function of sample translation in 
the transverse direction for a spot size of 2.2 //m, a 
bin size of 6 = 20, and three different NP 
concentrations with stretched exponential fits. As 
the concentration increases the width increases. 


concentrations of 1 wt%, 10 wt%, and 25 wt%, a 
bin size of 5 = 20, and a spot size of 2cr = 2.2 
jiim. We first measure the effect of the different NP 
concentrations on the coupling’s transverse transla¬ 
tion stability with the correlation coefficient data 
and stretched exponential fits shown in Figure [T3 
From Figure [T7| we find that the correlation curve 
width increases with increasing concentration, albeit 
only to a small degree. Table [Till lists the HWHMs 
for the three concentrations, from which we see that 
the increase in HWHM is larger than experimental 
uncertainty, but is relatively small compared to the 
changes due to the bin size and spot size. The in¬ 
crease in the translational stability with increasing 
concentration is due to the high concentration sam¬ 
ples having more uniform surfaces, such that there is 
less variation in the reflection matrix as the sample 
is translated in the transverse direction. 

While the effect of NP concentration is found to be 
small on the wavefront-sample coupling’s transverse 
translational stability, its influence is more appar- 
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Scaled Translation (units of z„) 

0.0 0.2 0.4 0.6 0.8 1.0 



FIG. 18: (Color online) Measured correlation 
coefficient as a function of sample translation along 
the optical axis for a spot size of 2.2 /im, a bin size 
of 6 = 20, and three different concentrations with 
stretched exponential fits. 


ent when we consider translation along the optical 
axis. Figure [18] shows the correlation coefficent as 
a function of 2 : position and stretched exponential 
fits for the three concentrations. From Figure [18] we 
find that the three curves intersect near the Rayleigh 
range (also observed for varying the bin size and 
spot size) and that as the concentration increases the 
HWHM (listed in Table [m]) of the curves decreases, 
implying that the lower concentration samples are 
more stable to translation along the optical axis. 

The mechanism behing the increased stability of 
low concentration samples to translation along the 
optical axis is related to the relative light penetra¬ 
tion depths of the different samples. As the con¬ 
centration increases, the scattering length decreases 
meaning that light penetrates less deeply into the 
sample. This depth can be viewed as defining an ef¬ 
fective interaction volume. For lower concentration 
samples this volume is larger than that of higher con¬ 
centrations. Therefore when translating the sample 
along the optical axis, the fractional change in the 
scattering volume will be smaller for large scatter¬ 
ing volumes than for small scattering volumes. This 
difference implies that the wavefront-sample inter¬ 
action will be more affected in the higher concentra¬ 
tion samples, than in the low concentration samples, 
which is observed. 

Lastly, we measure the effect of the NP concentra¬ 


tion on the wavefront-sample coupling’s rotational 
stability. Figure [T9] shows the correlation coefficient 
as a function of sample rotation for the three dif¬ 
ferent concentrations and fits to stretched exponen¬ 
tials. Similar to the transverse translation we find 
that higher concentrations result in increased rota¬ 
tional stability, with the HWHM, tabulated in Table 
mil increasing as a function of concentration. The ef¬ 
fect of increasing concentration is also found to be 
small, with the HWHM only varying by ~ 60 /irad 
between the 1 wt% and 25 wt% sample. This differ¬ 
ence is due to the surface being more uniform for the 
high concentration samples than for the low concen¬ 
tration samples. 

Apart from the correlation curve’s HWHM we also 
note from Figure [19] that the correlation offset (cq 
in Equation [12]) is found to increase with increasing 
concentration. This result is due to the nature of the 
speckle pattern as a function of concentration. As 
the concentration increases the speckle grain is found 
to decrease (i.e. individual speckles become smaller) 
resulting in a larger speckle density. If we compare 
two different regions of the speckle pattern we find 
that the similarity of the two regions increases with 
the speckle density. This implies that comparing 
two images with high speckle density will result in 
more similarity and therefore a larger background 
correlation coefficient. 



Angle (iirad) 


FIG. 19: (Golor online) Measured correlation 
coefficient as a function of sample rotation for a 
spot size of 2.2 /im, a bin size of 5 = 20, and three 
different sample concentrations with Gaussian fits. 
The half-width angle decreases as a function of 
increasing concentration, suggesting that higher 
concentrations are more stable under rotation. 
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Concentration 

X Half-Width (nm) 

Z Half-Width (nm) 

Rot. Half-Width (/xrad) 

1 

628 ± 17 

2582 ± 80 

215.6 ± 1.6 

10 

670 ± 22 

2090 ± 51 

237.1 ± 1.4 

25 

715 ± 16 

1365 ± 25 

276.4 ± 2.7 


TABLE III: Translational and rotational HWHMs for different NP concentrations 


V. CONCLUSIONS 

In the above sections we have explored theoreti¬ 
cally and experimentally the effects of three differ¬ 
ent system parameters (bin size, spot size, and NP 
concentration) on the wavefront-sample coupling’s 
translational and rotational stability for an opti¬ 
mized reflection system. We find that the simula¬ 
tions - using a model developed from the random 
phase Gaussian beam model xTa and Vellekoop’s 
random unitary matrix model - are within qual¬ 
itative agreement with experimental results. The re¬ 
sults of experiment and simulation are summarized 
below. 

When changing bin sizes it is found that larger 
bins provide more translational stability in both the 
transverse and z directions. This improvement is 
due to the fractional translation change being in¬ 
versely proportional to the bin size, which means 
that larger bins result in the on sample wavefront 
changing less during translation than for smaller 
bins. While the bin size affects the translational 
stability it does not affect the wavefront-sample cou¬ 
pling’s rotational stability. This is expected as: I) 
the equations describing the effect of sample rota¬ 
tion (Equations [2Ql[23|) are independent of the bin 
size and 2) for small angles the bin size is invariant 
under rotation. 

The next system parameter varied, the on sample 
spot size, is found to influence both the transverse 
translation and rotational stability of the wavefront- 
sample coupling, but not the stability of translation 
along the optical axis. As the spot size is increased 
the transverse translational stability is found to in¬ 
crease, while the rotational stability decreases. The 
benefit of larger spot sizes on the transverse stabil¬ 
ity is due to the fractional wavefront change (during 
transverse translation) being inversely proportional 
to the spot size. As the spot size increases the ef¬ 
fect of translating in the transverse direction has a 


smaller effect on the wavefront seen by the sample. 
However, in the case of rotation, the increased spot 
size results in larger wavefront changes due to Equa¬ 
tion [23Is dependence on the spot size. Einally, the 
lack of an effect on translation along the optical axis 
from changing the spot size is due to the 2 : transla¬ 
tional stability depending on the fractional bin size, 
and not the absolute bin size. While the absolute bin 
size changes with spot size, the fractional bin size is 
unchanged, thus explaining the observed effect. 

The last system parameter varied, the sam¬ 
ple’s NP concentration, is found to increase the 
wavefront-sample’s coupling stability to transverse 
translation and rotation when increased, while also 
decreasing the stability to translations along the op¬ 
tical axis. These results are due to higher concen¬ 
tration samples having more surface uniformity and 
smaller scattering volumes. 

The sum total of these results suggest that in prac¬ 
tice a balance needs to be struck in order to obtain 
the best possible positional stability. This obser¬ 
vation is due to both spot size and sample NP con¬ 
centration being found to produce conflicting results 
(e.g. smaller spot sizes increase rotational stabil¬ 
ity while decreasing transverse stability). Addition¬ 
ally, changing these parameters for stability needs 
to be weighed against their effects on the maximum 
possible enhancement, described in Ref. Eor 

instance, larger bin sizes are more positionally sta¬ 
ble but produce smaller enhancements. Therefore 
in some applications it may be more important to 
sacrifice stability for larger enhancements. 
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